In [7]:
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import scipy.sparse.linalg as sp
import sympy as sym
from scipy.linalg import toeplitz
import ipywidgets as widgets
import matplotlib as mpl
mpl.rcParams['font.size'] = 14
mpl.rcParams['axes.labelsize'] = 20
mpl.rcParams['xtick.labelsize'] = 14
mpl.rcParams['ytick.labelsize'] = 14
sym.init_printing()
Defining symbolic variables
In [9]:
xj, h, x = sym.symbols('xj h x', reals=True)
u1, u2, u3, u4, u5 = sym.symbols('u1 u2 u3 u4 u5', reals=True)
In [10]:
# The little Lagrange polynomials
def l_Lagrange(X,x,i):
l=1
n=len(X)
for k in np.arange(i):
l*=(x-X[k])
for k in np.arange(i+1,n):
l*=(x-X[k])
return l
# The Lagrange polynomials
def L_Lagrange(X,x,i):
num = l_Lagrange(X,x,i)
den = l_Lagrange(X,X[i],i)
L=num/den
return L
# The Lagrange interpolation
# X : x - data points
# x : the symbolic variable
# Y : y - data points
def P_Lagrange(X,x,Y):
P=0
n=len(X)
for i in np.arange(n):
P+=Y[i]*L_Lagrange(X,x,i)
return P
In [6]:
# x - data points
X=(0,1,2,3,4,5,6)
# y - data points
Y=(0,1,-2,1,0,2,7)
# The interpolation
P=P_Lagrange(X,x,Y)
# Just showing the polynomial. Do you see the Lagrange polynomials?
P
Out[6]:
In [11]:
# And a simplified version of the polinomial!
sym.simplify(P)
Out[11]:
Plotting what we got
In [12]:
# We need to plot in a finer grid
xx=np.linspace(min(X),max(X),1000)
# Here we create a lambda function for our interpolation
Pl=sym.lambdify(x,P)
# The next line vectorize the lambda function, it is not need here but it can be used
# Plv=np.vectorize(Pl)
plt.figure()
# Plotting the interpolated data
plt.plot(xx,Pl(xx),label='$P(x)$')
# Plotting the data points
plt.plot(np.array(X),np.array(Y),'k.',markersize=20,label='Data points')
plt.grid(True)
plt.ylim([np.min(Pl(xx))-1,np.max(Pl(xx))+1])
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend(loc='best')
plt.show()
In [21]:
# x - data points. Notice that they are all symbolic variables!
X=(xj-2*h,xj-h,xj,xj+h,xj+2*h)
# y - data points. Also symbolic variables!
Y=(u1,u2,u3,u4,u5)
## In-class
#X=(xj-h,xj,xj+h)
#Y=(u1,u2,u3)
# The Lagrange interpolation still works!
P=P_Lagrange(X,x,Y)
P
Out[21]:
In [22]:
# Here we compute the derivative of the polynomian P(x), i.e. we obtain P'(x) and we definy it as Pp (P-prime).
Pp=sym.diff(P,x)
Pp
Out[22]:
In [23]:
# Now we evaluate P'(x) at x=x_j.
PpE=Pp.subs({x:xj})
PpE
Out[23]:
In [ ]:
In [ ]: